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MODELLING  THE  RESPONSE  OF  THICK  COMPOSITE  MATERIALS 
DUE  TO  AXISYMMETRIC  SHOCK  LOADING 


I.  INTRODUCTION 

The  potential  use  of  thick  composite  materials  for  future  naval  vessels  requires 
an  understanding  of  their  response  to  dynamic  loadings  such  as  those  arising 
from  impact  and  underwater  explosives.  Although  numerous  investigations 
have  been  conducted  on  the  linear  vibration  of  composite  plates  [1],  little  has 
been  done  in  the  area  of  response  to  extreme  transient  loadings,  particularly 
when  nonlinear  behavior  associated  with  material  damage  is  involved.  Some 
work  in  the  area  of  material  response  to  uniaxial-strain  shock  loading  has  been 
done  [2-6],  including  identification  of  spallation  thresholds.  Response  of  the 
thick  composite  due  to  multi-dimensional  deformation  type  shock  loading  has 
received  much  less  attention. 

In  this  report,  we  report  on  our  preliminary  efforts  to  develop  a  methodology 
for  predicting  the  response  of  thick  composite  materials  to  multidimensional 
shock  loadings.  Initially  we  are  limiting  our  consideration  to  composite 
material  layups  that  result  in  transversely  isotropic  behavior  in  the  plane  of  the 
fibers.  This  assumption  along  with  consideration  of  loadings  that  are 
symmetric  about  the  axis  of  material  symmetry  greatly  simplifies  the  fully 
three-dimensional  problem,  while  still  being  of  both  scientific  and  practical 
interest.  The  configuration  described  above  is  depicted  schematically  in 
Figure  1.1.  A  principal  practical  interest  in  this  configuration  arises  in  the 
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study  of  the  explosive  bulge  test  (EBT).  The  EBT  has  been  used  over  a  long 
history  [7]  to  study  the  deformation  and  fracture  behavior  of  ship  steels  under 
high  rate  loadings,  particularly  the  behavior  of  weldments  [8].  The  use  of  the 
test  continues  today  [9-10].  Several  current  studies  are  underway  to  use  a  form 
of  the  EBT  to  study  the  response  of  composite  materials  to  shock  loading 
[11,12]. 

Although  a  complete  description  of  the  EBT  is  not  part  of  the  s^^ope  of  this 
paper,  it  will  briefly  be  described  here.  The  majority  of  the  tests  are  conducted 
underwater  since  the  water  can  effectively  transmit  the  shock  to  the  plate 
surface,  and  because  of  its  apparent  similarity  to  actual  conditions  experienced 
by  structural  components  subjected  to  underwater  explosives.  The  plate  to  be 
tested  is  supported  on  an  anvil  having  a  circular  hole  in  a  manner  such  that  the 
rear  surface  of  the  plate  is  evacuated.  An  explosive  is  placed  at  a  given 
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standoff  from  the  submerged  side  of  the  plate  and  detonated  electrically. 
Instrumentation  used  during  the  tests  varies  somewhat,  but  recording  of 
deflection  histories  of  the  plate,  pressure  histories  at  the  front  surface,  and 
strain  histories  at  the  rear  surface  are  common. 

The  material  and  loading  symmetries  considered  in  this  paper  have  practical 
interest  beyond  the  consideration  of  the  EBT.  In  cases  where  actual 
configuration  of  naval  components  is  such  that  the  layup  is  transversely 
isotropic  (or  nearly  so),  the  local  response  to  such  events  as  torpedo  impact  or 
shock  wave  due  to  a  nearby  explosive  can  be  appropriately  treated  using  the 
configuration  considered  here  as  long  as  the  offset  distance  is  much  smaller 
than  the  radius  of  the  shell  being  considered  and  the  events  of  interest  occur 
during  times  such  that  the  outer  boundaries  can  be  appropriately  considered. 

In  Section  H,  the  kinematics  associated  with  the  configuration  discussed 
above  are  developed  along  with  the  equations  of  motion  in  the  cylindrical 
coordinate  system  specialized  for  axisymmetric  deformation.  Section  HI 
contains  the  development  of  the  material  constitutive  model.  In  this 
preliminary  investigation,  the  laminated  composite  material  is  idealized  as  a 
transversely  isotropic,  homogeneous  linear  solid.  The  subsequent 
development  of  a  nonlinear  damage  model  and  its  applications  are  presented 
in  later  sections.  Material  properties  are  discussed  in  Section  IV.  The 
computations  presented  here  have  been  performed  using  graphite/epoxy  as  the 
model  material,  both  because  of  present  interest  in  the  material  and  its  fairly 
extensive  characterization  at  least  under  static  loading  conditions.  Under  high 
loading  rates,  however,  experimental  data  on  the  material  behavior  is  sparse 
[13,14].  Section  V  provides  a  brief  discussion  of  loading  associated  with 
underwater  explosives.  For  the  purposes  of  this  report,  use  is  made  of 
empirical  formulas  to  determine  amplitude  and  shape  of  the  shock  wave 
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loading.  A  discussion  of  the  numerical  method  used  is  included  in  Section  VI. 
Explicit  integration  of  the  balance  equations  of  mass  and  momentum  is 
performed  using  the  computer  code  PRONTO-2D  [15].  Several  sample 
problems  are  considered  in  Section  Vn  to  illustrate  the  methodology 
considered  for  shock  response  of  thick  composites.  The  development  of  the 
damage  model  is  contained  in  Section  VUI  and  several  cases  of  homogeneous 
deformation  are  considered  in  Section  DC.  Finally,  preliminary  application  of 
the  damage  model  to  the  multi-dimensional  response  to  shock  loading  is 
considered  in  Section  X. 
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n.  KINEMATICS  AND  EQUATIONS  OF  MOTION 


The  kinematics  to  describe  the  defonnation  associated  with  the  configuration 
shown  in  Figure  1.1  is  presented  in  this  section.  The  reference  configuration 
of  the  body  is  defined  by  particles  that  occupy  position  X  in  that 
configuration.  The  current  configuration,  x,  is  defined  by  the  motion,  such 
that 


x  =  X(X,t).  (2.1) 

Taking  the  time  derivative  of  (2.1),  holding  X  fixed,  defines  the  material 
velocity 

y  =  j  =  (^X  (X,  t) )  .  (2.2) 

The  velocity  gradient 


L  =  ^ 

^  ax 

can  be  decomposed  into  symmetric  and  anti-symmetric  parts 


(2.3) 


L  =  P  +  W, 


(2.4) 


where  the  skew  symmetric  tensor  W  is  the  spin  tensor.  For  small  strains  p  -  e 


so  that  the  strain  tensor  is  defined  by 


t 

e  =  jpdt’ .  (2.5) 

0 

Referring  to  cylindrical  coordinates,  X  =  (X’,  X^x’)  =  (Z,  R,  0)  and 
x=  (x^x^x^)  =  (z,r,  0),  for  axisymmetric  deformation  where  the  motion  is 
independent  of  0,  %  =  X(Z, R,t)  .  The  strain  tensor,  therefore  has  non-zero 
components  Eu,  £22.  £33.  and  e|2  =  £21.  all  of  which  are  independent  of  0. 


The  balance  of  mass  requires  that 


(2.6) 


where  p(x,t)  and  Pr(X)  are  the  mass  densities  in  the  current  and  reference 
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configuration  respectively  and  J  is  the  Jacobian  of  the  deformation.  Balance 
of  linear  momentum  leads  to  the  equation  of  motion 


diva  +  pb  =  p^^. 


(2.7) 


which,  for  axisymmetric  deformation  (neglecting  the  body  force  b),  reduces 


.  ^rz  .  1 


3f  ■^7<°rr-®6e>  -  Pdt 


^^rz^^^zz^l  _  I'-z 

3r  3z  r*^rz  ^dt 


(2.8) 
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in.  CONSTITUTIVE  MODEL 


Linear  Elastic  Model 

The  generalized  Hooke’s  Law  relating  stress  and  strain  can  be  written  as 


where  C  is  the  fourth-order  elasticity  tensor,  which  for  completely  anisotropic 
materials  contains  21  independent  components.  Materials  containing  three 
planes  of  symmetry  are  termed  orthotropic  and  contain  9  independent 
components.  In  addition,  if  there  is  one  plane  in  which  mechanical  properties 
are  equal  in  all  directions,  the  material  is  transversely  isotropic  and  the  number 
of  independent  material  constants  is  reduced  to  five. 


The  composite  material  under  consideration  is  idealized  as  transversely 
isotropic  thus  requiring  specification  of  only  the  five  components. 
Considering  the  axisymmetric  deformation  discussed  in  Section  II,  Equation 
(3.1)  can  be  written  as 


^11 

^11^12^12  ® 

^11 

^22 

^12  ^22  ^23  ® 

^22 

^33 

C12  C23  C22  0 

^33 

512 

0  0  0  C44 

.^12 

(3.2) 


where  the  contracted  notation  has  been  used  for  the  fourth  order  tensor  C.  The 
components  of  C  given  in  terms  of  the  engineering  properties  are 

C  =  (1  -V23)Eii 

"  >-''23-(2v?2E22)/E„ 

”  (l+Vjj)  [1-V23-(2v2jE22)/E„] 


c 


12  " 


12^22 

1  -  V23  -  (2v  j2E22^  1 


(3.3) 
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c  -  ^^23  ^^12^22)  ^^11^^22 

(l+Vn-V23-(2v22E22)/E„] 

C44  =  2Gj2  . 

For  modelling  isotropic  materials,  only  two  independent  material  properties 
are  required.  Thus,  specification  of  the  Young’s  modulus,  E,  and  the  Poisson’s 
ratio,  V,  is  sufficient  to  characterize  the  linear  elastic  behavior  of  the  material. 


G  = 


E 

2(H-v) ' 


(3.6) 


For  studying  problenis  involving  fluid-structure  interaction,  modelling  the 
fluid  using  finite  elements  offers  many  advantages  over  using  boundary 
element  techniques,  such  as  doubly  asymptotic  approximations  (DA A)  [17]. 
Material  nonlinearity  as  well  as  cavitation  can  be  addressed  fairly  easily  using 
finite  elements,  while  either  can  b".  modelled  only  with  great  difficulty  or  not 
at  all  with  boundary  element  techniques.  In  studies  of  material  damage  it  is 
expected  that  pressures  acting  in  the  fluid  will  be  sufficiently  high  that  the 
nonlinear  fluid  behavior  is  significant.  In  addition,  shock  waves  propagating 
in  a  fluid  and  impinging  on  a  flat  plate  in  the  configuration  of  Figure  1.1  can 
be  expected  to  result  in  cavitation  as  waves  reflect  back  from  the  free-rear 
surface  of  the  solid.  The  use  of  non-reflecting  boundaries  [18,19],  discussed 
in  more  detail  in  Section  VI,  limits  the  amount  of  fluid  that  must  be  modelled, 
which  is  often  cited  as  a  drawback  in  using  finite-elements  to  model  fluid 
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behavior. 


The  constitutive  model  for  the  fluid  is  taken  as  a  simple  hydrodynamical 
model  where 

Ojj  =  -p6. j  =  Kgri  ( 1  +  K jTi  +  ¥.2^  +  . . . )  for  (T|  >  0)  (3.7) 

where  Tj  =  1-  po/p  and  p  =0  for  ti  <0. 
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rV.  MATERIAL  PROPERTIES 


For  the  simulations  conducted  in  this  report  the  model  material  used  is  for  a 
graphite/epoxy  with  a  transversely  isotropic  layup.  Tsai  [16]  gives  properties 
for  a  graphite/epoxy  composite  composed  of  T300  fiber  with  N5208  matrix. 
The  properties  given  are  for  p,  E,  and  v,  which  are  the  in-plane  averaged 
properties  corresponding  to  E22  and  V23,  in  the  notation  given  in  Section  HI. 
The  properties  given  are 

p  =  1600  kg/m^ 

E22  =  69.0x  lO^MPa  (4.1) 

V23  —  0.30 . 

The  through-thickness  modulus  E22,  the  transverse  shear  modulus,  G22,  and 
the  V22  Poisson’s  ratio  are  assumed  to  be  equal  to  the  transverse  modulus,  the 
in-plane  shear  modulus,  and  the  Poisson’s  ratio,  respectively,  of  a  single  ply 
given  in  [16],  such  that 

Ell  =  10.3x10^  MPa 

Gi2  =  7.17x10^  MPa  (4.2) 

V21  =  0.28 . 

For  comparative  simulations  involving  steel  the  properties  used  are 

p  =  7779  kg/m^ 

E  =  206.8  X  10^  MPa  (4.3) 

V  =  0.30 . 

Material  properties  for  water,  referring  to  the  constitutive  model  of  Section  HI, 
are  given  by  (with  K2=0) 

p  =  998  kg/m^ 

Ko  =  2.33x  lO^MPa  (4.4) 

Ki  =  1.75. 
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V.  LOADING  CONDITIONS  CONSIDERED 


The  theory  of  shock  waves  from  underwater  explosion  has  been  extensively 
treated  and  is  defined  primarily  by  the  Kirkwood-Bethe  theory  [18].  Because 
of  its  complexity,  however,  Kirkwood-Bethe  is  infrequently  used,  and  it  is 
replaced  for  the  most  part  by  semi-empirical  scaling  laws  such  as  the  Arons 
[19],  which  gives  the  peak  pressure  of  the  shock  wave  in  terms  of  weight  of 
TNT,  W,  and  distance  from  the  center  of  charge,  R.  Pressures  from  other  types 
of  explosive  are  commonly  determined  by  defining  them  in  terms  of  an 
equivalent  weight  of  TNT.  The  scaling  law,  which  has  been  verified 
experimentally  over  several  orders  of  magnitude  is  given  by 


Pm  =  Kp  (w'-^/R) 


(5.1) 


When  R  is  given  in  meters  and  W  in  kilograms,  Kp  =  5.19x10^  and  is  given 
in  kilopascals.  When  R  is  given  in  feet  and  W  in  pounds,  Kp  =  2.16x10'^  and 
Pni  is  given  in  pounds  per  square  inch.  The  shape  of  the  shock  wave  is  taken 
to  have  an  instantaneous  rise  with  an  exponential  decay.  The  pressure  is  then 
given  by 


-J. 

P  (t)  =  P„,e  *'  (5.2) 

where  arrival  of  the  wave  is  taken  to  be  at  t  =  0.  The  characteristic  time,  Sj,  is 
also  given  by  an  Arons  scaling  law 


e  =  K  W^  W^/R 

I  ^ 


(5.3) 


where  Gj  is  given  in  microseconds.  K^  =  92.5  when  R  is  given  in  meters  and  W 
in  kilograms  and  K^  =  58  when  R  is  in  feet  and  W  in  pounds.  Rogers  [20] 
reports  that  (5.3)  has  been  experimentally  verified  over  a  smaller  range  of 
values  of  (W^'^  /R).  Problems  in  using  the  empirical  forms  given  by  (5.1)  and 
(5.3)  are  also  noted  [10]  for  small  offsets.  Nevertheless  use  of  the  empirical 
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forms  are  sufficient  to  illustrate  the  modelling  considered  in  this  report, 
although  caution  should  be  used  if  these  empirical  forms  are  used  in  trying  to 
correlate  computed  results  with  experimental  data.  A  range  of  explosive 
weights  and  offset  distances  is  considered  in  giving  the  peak  pressures  and 
characteristic  times  in  Tables  5.1  and  5.2,  developed  from  (5.1)  and  (5.3).  As 

TABLE  5.1  PEAK  PRESSURE  (PSI)  FOR  DIFFERENT  EXPLOSIVE 
WEIGHTS  AND  OFFSET  DISTANCE 

R  W  (LBS) 


(FT) 

in 

5n 

100 

5QQ 

IDQQ 

2000 

3000 

2 

23494 

43076 

55927 

- 

- 

- 

- 

5 

8342 

15295 

19859 

36411 

47273 

61377 

71504 

10 

3812 

6989 

9074 

16637 

21600 

28044 

32672 

20 

\1A1 

3193 

4146 

7601 

9869 

12814 

14928 

50 

618 

1134 

1472 

2699 

3504 

4550 

5301 

100 

283 

518 

673 

1233 

1601 

2079 

2422 

200 

129 

231 

307 

564 

732 

950 

1107 

seen  from  the  tables  the  peak  pressures  vary  tremendously  but  the 
characteristic  time  is  almost  entirely  in  the  range  of  lO"^  to  10'^  seconds. 
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TABLE  5.2  CHARACTERISTIC  TIME  0t  (MSEC)  FOR  DIFFERENT  EXPLOSIVE 

WEIGHTS  AND  OFFSET  DISTANCE 


R 

(FT) 

W(LBS) 
too  5QQ 

1000 

2000 

3000 

2 

.1229 

.1868 

.2237 

- 

- 

- 

- 

5 

.1504 

.2285 

.2737 

.4159 

.4980 

.5963 

.6626 

10 

.1752 

.2662 

.3187 

.4844 

.5800 

.6945 

.7718 

20 

.2040 

.3100 

.3712 

.5641 

.6756 

.8089 

.8989 

50 

.2496 

.3793 

.4542 

.6901 

.8264 

.9896 

1.0996 

100 

.2907 

Mil 

.5290 

.8038 

.9626 

1.1526 

1.2808 

200 

.3386 

.5145 

.6161 

.9362 

1.1211 

1.3425 

1.4918 
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VI.  NUMERICAL  METHODS 


Solution  of  the  balance  equations  of  mass  and  momentum  for  the 
configuration  shown  in  Figure  1.1  is  performed  numerically,  using  the  finite- 
element  computer  code  PRONTO-2D  [15].  PRONTO-2D  is  a  two- 
dimensional,  explicit  integration,  finite-element  cr.de,  developed  specifically 
for  solution  of  transient  problems.  The  code  is  developed  in  a  Lagrangian 
formulation  [21].  Simple  four-noded  quadrilateral  uniform  strain  elements, 
with  hourglass  control  [22]  are  used  for  the  discretization. 

In  computing  the  response  to  shock  loading,  that  is,  a  loading  that  produces  a 
discontinuity  in  the  particle  velocity,  artificial  viscosity  [23]  must  be  applied 
to  the  numerical  procedure.  The  effect  of  artificial  viscosity  is  to  smear  a 
shock  front  across  several  elements,  thus  replacing  the  shock  front  with  one 
that  has  a  rapid  but  finite  rise  time.  The  effect  of  this  smearing  of  the  shock  on 
the  solution  is  discussed  further  in  the  following  sections. 

The  numerical  procedure  employed  also  makes  use  of  nonreflecting  boundary 
conditions  [24,  25].  These  boundaries  are  particularly  useful  for  limiting  the 
extent  of  fluid  to  be  modelled.  Use  of  a  simple  free  surface  would  require  that 
a  sufficient  volume  of  fluid  be  modelled  such  that  reflections  off  that  surface 
do  not  reflect  and  impinge  again  on  the  body  of  interest.  The  use  of 
nonreflecting  boundaries,  however,  are  intended  to  behave  as  if  the  fluid 
extended  to  infinity  while  limiting  the  actual  region  to  be  modelled  [25]. 
Implementation  of  the  nonreflecting  boundary  for  the  fluid  is  accomplished  by 
applying  tractions  normal  to  the  boundary  such  that  the  stresses  at  the  surface 
cancel.  The  normal  stress  applied,  a„,  is  determined  from 

<J„  =  pCpV„  (6.1) 

where  p  is  the  current  fluid  density,  Cp  is  the  bulk  sound  speed  in  the  fluid  and 
Vj,  is  the  fluid  particle  velocity  in  the  direction  normal  to  the  boundary. 
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The  fluid-solid  interface  is  modelled  using  a  contact  surface  algorithm  in 
which  both  media  are  deformable.  The  algorithm  uses  kinematic  constraints, 
in  which  the  accelerations  of  the  nodes  along  the  surface  are  modified.  Final 
accelerations  are  computed  such  that  one  surface  may  not  penetrate  the  other, 
while  at  the  same  time  establishing  constraint  forces  that  conserve 
momentum.  At  each  time  step  the  motion  is  first  computed  without  the 
kinematic  constraints  to  determine  the  extent  of  contact  between  the  surfaces. 
From  this  computed  motion,  the  penetration  forces  necessary  to  satisfy  the 
kinematic  constraints  are  calculated  and  used  to  modify  the  accelerations. 


15 


Vn.  NUMERICAL  RESULTS  FOR  SHOCK  LOADING 


In  this  section  several  sample  problems  are  considered  in  order  to  illustrate  the 
methodology  considered  for  addressing  the  shock  response  of  thick 
composites  under  multi-axial  loading.  The  model  composite  material  used  is 
graphite/epoxy  with  the  properties  given  in  Section  IV.  Several  comparisons 
with  steel  are  presented.  Where  possible  comparisons  are  also  made  with 
analytical  solutions.  The  problems  are  intended  to  illustrate  various  phases  of 
the  methodology  considered  here  for  evaluating  multi-dimensional  material 
response  to  shock  loading,  such  as  fluid-structure  interaction,  wave 
propagation  through  the  material,  development  of  tensile  stresses  in  the 
through-thickness  direction  of  the  composite  and  finite-element 
discretization. 


.Fdan!LWay£ 

The  first  problem  considered  is  the  response  of  a  25.4  mm  thick  infinite  plate 
subjected  to  a  plane  wave  of  temporal  distribution  given  by  Equation  (5.2). 
Since  we  are  interested  in  normalizing  the  response  only  the  linear  behavior 
of  the  fluid  is  considered,  i.e  Ki=K2  =  0  in  Equation  (3.7).  The  characteristic 
time,  6t,  is  taken  to  be  0.0001  sec.  The  finite  element  model  for  the  problem  is 
shown  in  Figure  7.1.  The  model  uses  a  mesh  of  19x19  elements  for  the  fluid, 
which  is  254  mm  in  the  z-direction,  and  39x19  elements  for  the  plate.  The 
mesh  is  considered  to  be  fairly  coaree  and  is  used  here  for  illustration 
purposes.  The  sides  given  by  r=0  and  r=ri  are  constrained  such  that  the 
velocity  v^  =0.  The  plate  is  free  in  the  z-direction.  A  uniform  spatial  pressure 
is  applied  to  the  outer  boundary  surface  of  the  fluid  and  the  silent  boundary 
condition  is  also  applied  to  this  boundary  to  prevent  reflection  of  the  waves 
back  from  the  boundary  surface  after  it  has  reflected  off  the  plate.  These 
boundary  and  loading  conditions  result  in  a  uniaxial  deformation  field,  which 
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could  more  appropriately  be  described  using  a  one-dimensional  finite- 
difference  technique.  The  purpose  here,  however,  is  to  evaluate  the  use  of  the 
two-dimensional  finite  element  modelling  technique  by  comparing  results 
with  an  analytic  solution.  The  uniaxial  deformation  allows  the  problem  to  be 
modelled  independent  of  the  size  of  the  mesh  in  the  r-direction,  thus  a  fairly 
narrow  slice  is  used  as  shown  in  Figure  7.1. 


ilic  analytical  solution  for  the  problem  described  above  can  be  obtained 
directly  from  consideration  of  transient  elastic  plane  waves  in  an  elastic 
media.  The  wave  profile  described  by  Equation  (5.2)  travels  undisturbed 
through  the  fluid  at  a  velocity  equal  to  its  sound  speed  until  it  encounters  the 
boundary  of  the  plate.  At  that  point  a  portion  of  the  wave  is  transmitted  and  a 
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portion  of  it  is  reflected.  The  relative  fractions  of  each  is  governed  by  the 
density  and  wave  velocity  of  the  two  media.  At  the  boundary  between  the  two 
media  the  normal  stresses  must  be  equal  at  every  instant  of  time  and  the 
normal  particle  velocities  must  be  equal.  This  holds  in  compression,  but  when 
tension  develops  across  the  interface,  separation  of  the  fluid  from  the  solid 
occurs.  The  boundary  conditions  result  in  the  following  relationships  between 
the  transmitted,  reflected  and  incident  stress  wave  magnitudes  (Oj,  Or,  and  O], 
respectively) 

0^.=  [(2pV)/(pV  +  pc)]Oj  (7.1a) 

=  [(pV-pc)/(pV  +  pc)]Oj  (7.1b) 

where  p  is  the  density  and  c  is  the  wave  speed.  The  primes  denote  the  second 
medium,  which  in  this  case  is  the  solid,  and  the  unprimed  quantities  are  for  the 
fluid.  The  wave  speed  for  the  fluid  is  given  by 

c  =  (Ko/p)>''2.  (7.2) 

The  wave  velocity  for  a  transversely  isotropic  solid  is  given  by 

c=(C„/p)>/2  (7.3) 

or  in  the  case  of  an  isotropic  solid  by 

c  =  [(>.  +  2G)/p]^/2  (7  4) 

The  stress  in  the  fluid  is  then  the  sum  of  the  incident  and  reflected  waves  and 
the  stress  in  the  solid  is  that  given  by  the  transmitted  wave.  Examination  of 
Equation  (7.1a),  in  view  of  the  material  properties  given  in  Section  IV,  shows 
that  1 .93  times  the  incident  wave  would  be  transmitted  to  steel,  but  only  1 .46 
times  the  incident  wave  would  be  transmitted  to  the  graphite/epoxy. 
Implications  of  this  on  the  ability  of  the  material  to  withstand  shock  loading 
are  discussed  later. 

The  transmitted  wave  propagates  through  the  solid  until  it  encounters  the  free 
rear  surface.  Conditions  at  this  boundary  can  be  obtained  from  Equations 
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(7.1a)  and  (7.1b)  by  considering  the  first  medium  to  be  the  solid  (unprimed) 
and  the  second  medium  to  have  p'c'  =0.  Thus  from  Equations  (7.1a)  and 
(7.1b),  the  transmitted  wave  is  equal  to  zero  and  the  reflected  wave  is  equal  to 
-  Oj.  This  wave  is  superimposed  with  the  incident  wave  in  the  solid  to  give  the 
total  stress  in  the  solid.  By  applying  these  boundary  conditions  for  each 
interface  encountered,  the  stress  history  in  both  the  solid  and  the  fluid  can  be 
obtained.  These  analytical  results  are  shown  and  compared  with  the  finite 
element  calculations  to  illustrate  the  effects  of  discretization  and  artificial 
viscosity  on  numerical  accuracy. 

Results  of  the  finite-element  analysis  are  given  in  terms  of  pressure  histories 
in  the  fluid  or  stress  histories  in  the  solid  at  specific  points  and  then  compared 
with  the  analytical  results.  The  first  simulation  is  for  the  plane  wave  impinging 
on  the  graphite/epoxy  plate.  Figure  7.2  shows  the  computed  pressure  history 
at  the  element  adjacent  to  the  centerline  (R=0)  and  furthest  out  in  the  fluid  (Z 
=  -247.00  mm).  The  idealized  shock  front  now  has  a  finite  rise  and  the  decay 
exhibits  some  oscillations,  which  are  numerical  in  origin.  The  rise  that  occurs 
after  0.0003  seconds  is  a  result  of  the  reflected  wave  from  the  front  of  the 
plate.  The  analytical  solution  for  the  same  point  in  space  is  also  given  in  Figure 
7.2.  The  computed  solution  agrees  well  with  the  analytic  solution  concerning 
the  magnitude  of  the  initial  wave  as  well  as  the  slope  of  the  decay.  The 
magnitude  of  the  reflected  wave  is  under  predicted  by  the  numerical  solution. 
Figure  7.3  shows  the  pressure  history  for  the  element  at  the  middepth  (Z  =  - 
1 13.(X)  mm)  of  the  fluid.  Spreading  of  the  shock  front  is  more  pronounced  as 
are  the  numerical  oscillations  in  the  decay.  The  rise  at  0.00025  seconds  is  not 
numerical,  but  is  due  to  the  reflection  from  the  face  of  the  plate.  The 
comparable  analytic  solution  is  also  given  in  Figure  7.3.  Comparison  between 
the  analytic  and  numerical  solutions  is  reasonable. 
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Figure  7.3  Finite<elenient  and  analytical  results  for  pressure  in  the  fluid 
at  z=-113  mm  due  to  plane  wave  shock  loading  on  a  graphite/epoxy 
plate. 
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The  history  of  the  normal  component  of  stress  in  the  z-direction  in  the  solid  at 
the  element  closest  to  the  interface  (Z  =  0.67  mm)  with  the  fluid  is  shown  in 
Figure  7.4.  The  analytic  solution  is  also  given  in  Figure  7.4.  Comparison  of 
the  two  shows  that  the  finite-element  solution  under  predicts  the  peak 
compressive  stress  reached  at  the  surface  of  the  plate.  The  duration  of  the 
compressive  pulse,  which  is  directly  related  to  the  travel  time  through  the 
thickness  of  the  plate,  is  longer  in  the  numerical  solution.  The  reason  for  both 
of  these  findings  can  be  traced  to  the  spreading  of  the  shock  rise  as  the  wave 
propagates  through  the  fluid  and  subsequently  the  solid.  The  similar  stress 
history  for  a  point  at  mid-thickness  of  the  plate  (Z  =  14.04  mm)  is  shown  in 
Figure  7.5  for  the  finite-element  and  analytic  solution,  respectively.  Here  the 
under  prediction  of  the  compressive  stress  in  the  solid  is  quite  significant  and 
also  significant  is  the  fact  that  the  analytical  solution  is  predicting  the 
existence  of  fairly  large  tensile  stresses  that  are  almost  negligible  in  the 
numerical  solution.  As  seen  in  comparing  the  two  results  is  that  the  frequency 
of  the  oscillations  is  quite  similar  between  the  two  with  the  principal 
difference  being  that  the  numerical  solution  is  rounding  off  the  wave  profiles 
that  are  sharp  and  distinct  in  the  analytical  solution. 

For  the  purposes  of  comparison,  a  second  simulation  conducted  using  steel  as 
the  plate  material  is  presented.  Pressure  history  results  in  the  fluid  from  both 
the  finite-element  analysis  and  the  analytical  solution  are  given  in  Figure  7.6 
and  Figure  7.7.  The  principal  difference  to  note  in  the  pressure  histories  at  the 
two  points  in  the  fluid  is  the  much  larger  reflected  wave  that  occurs  at  both 
points  in  the  fluid.  The  stress  history  near  the  plate  interface,  in  Figure  7.8 
shows  that  the  maximum  compressive  stress  reached  is  higher  than  that  in  the 
graphite/epoxy  as  discussed  above.  The  computed  stress  at  midplane,  shown 
in  Figure  7.9  significantly  under  predicts  the  maximum  compressive  stress 
shown  in  the  analytical  solution.  This  is  more  severe  than  in  the  case  of  the 
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Figure  7.4  Finite-element  and  analytical  results  for  near  the  solid- 
fluid  interface  for  the  graphite/epoxy  plate  subjected  to  plane  shock 
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Figure  7.8  Finite-element  and  analytical  results  for  Oxz  near  the  solid- 
fluid  interface  for  the  steel  plate  subjected  to  plane  shock  wave  loading. 
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graphite/epoxy  plates  due  to  the  higher  wave  speed  in  the  steel.  The  spread  out 
shock  front  that  occurs  in  the  fluid  travels  faster  through  the  steel  than  in  the 
graphite/epoxy.  Thus  the  reflected  wave  from  the  free  rear  surface  of  the  plate 
has  already  reached  the  midplane  before  the  spread  out  compressive  wave  has 
been  able  to  reach  its  peak  value  as  shown  in  Figure  7.9. 

Tlie  results  of  the  comparisons  show  that  the  finite-element  analysis  correctly 
predicts  the  trends  seen  in  the  analytical  solution.  However  because  of  the 
coarse  discretization  employed  in  the  modelling,  spreading  of  the  shock  wave 
occurs  which  results  in  under  predicted  stresses  in  the  plate,  which  in  some 
cases  is  quite  significant.  It  should  be  noted  here,  however,  that  the  analytic 
solution  is  based  on  a  linear,  elastic,  homogeneous  solid.  Therefore  its 
application  to  wave  propagation  in  heterogeneous  media  cannot  account  for 
dispersion  that  occurs  as  a  result  of  differences  in  the  material  behavior  of  the 
two  phases.  In  fact  the  dispersion  that  takes  place  in  the  real  composite 
material  is  somewhat  similar  to  the  numerical  dispersion  that  occurs  as  a  result 
of  the  artificial  viscosity  and  discretization  used  in  the  analysis  as  discussed  in 
Section  VI.  Because  of  this  similarity  a  viscosity  term  similar  to  artificial 
viscosity  is  sometimes  employed  to  mimic  the  dispersion  seen  in  these 
materials  [26].  It  appears  that  a  finer  discretization  of  the  fluid  will  be  required 
to  prevent  spreading  of  the  shock  front,  particularly  when  the  stresses  in  the 
through  thickness  direction  are  of  interest.  In  subsequent  computations  a 
discretization  of  the  fluid  is  used  having  a  dimension  in  the  z-direction  that  is 
one  fourth  that  used  here. 

Spherical  Wave  on  an  Infinite  Plate 

The  second  problem  considered  is  for  a  spherical  wave  impinging  on  an 
infinite  plate.  The  temporal  distribution  of  the  wave  is  that  given  by  Equation 
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(5.2).  The  simulation  is  conducted  for  steel  since  analytical  results  are 
available  for  only  isotropic  materials.  A  25.4  mm  thick  plate  is  considered  and 
the  center  of  the  spherical  wave  is  taken  as  R=0  and  2^-257  mm.  The 
analytical  solution  is  given  by  Huang  [27],  and  utilizes  the  Mindlin  [28]  plate 
equations  of  motion  including  rotatory  inertia  and  transverse  shear  terms.  The 
plate  theory  considers  the  normal  stress  in  the  z-direction  to  be  zero,  which 
differs  from  the  finite-element  modelling  using  the  two  dimensional  uniform 
strain  quadrilaterals,  which  also  permits  higher  order  shear-deformation 
through  the  overall  thickness  than  the  Mindlin  theory.  Nevertheless,  the 
comparison  of  the  finite-element  solution  to  that  given  by  the  analytic  solution 
is  valuable  to  determine  the  ability  of  the  modelling  procedure  to  predict  plate 
response,  which  is  predominately  flexural  behavior.  This  represents  the 
extreme  from  the  problem  considered  in  the  preceding  subsection,  in  which  no 
flexural  behavior  existed  and  the  response  was  dominated  by  through 
thickness  wave  propagation. 

The  finite-element  model  for  this  problem  is  shown  in  Figure  7.10.  Radially 


Figure  7.10  Finite-element  discretization  for  impingement  of  an 
infinite’  plate  by  a  spherical  shock  wave. 


the  problem  is  modelled  out  to  R  =  420  mm.  Thus  the  model  can  be  used  to 
represent  an  ‘infinite  plate’  only  until  such  time  that  disturbances  can 
propagate  to  the  outer  boundary  of  the  plate  and  back  to  the  point  of  interest. 
Using  the  wave  velocity  computed  from  Equation  (7.4)  for  steel,  this  results 
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in  a  time  of  0.1448  msec  for  the  center  of  the  plate  (R=0).  The  portion  of  the 
fluid  modelled  extends  to  Z  =  -63.5  mm  and  uses  a  mesh  of  19  x  59  elements. 
The  solid  uses  a  mesh  of  14  x  99  elements. 

Results  are  given  in  terms  of  radial  stress  at  the  back  surface  of  the  plate 
normalized  by  the  peak  pressure  of  the  incoming  spherical  wave.  The 
history  is  given  in  terms  of  normalized  time,  x  =ct/h,  where  c  is  the  wave  speed 
in  the  fluid,  t  is  time  after  the  wave  impinges  on  the  plate,  and  h  is  the  plate 
thickness.  The  finite-element  solution  is  valid  for  an  infinite  plate  up  to  x  = 
17.4.  The  computed  and  analytical  results  are  shown  in  Figure  7.11.  The 


Ql - 1 - 1 - 1 

0  5  10  15 

ct/h 


Figure  7.11  Computed  and  analytical  results  for  radial  stress  in  an 
infinite  plate  subjected  to  a  spherical  shock  wave. 

comparison  shows  good  overall  agreement,  particularly  with  respect  to  the 
peak  radial  stress  reached.  Discrepancies  are  expected  since  the  finite-element 
solution  contains  through-thickness  stress  and  the  analytical  solution  does  not. 


Simply  Supported  Plate 


In  tWs  subsection  the  response  of  a  simply  supported  plate  to  a  spherical  shock 
wave  is  considered.  Plates  of  25.4  mm  and  50.8  mm  thicknesses  using 
properties  for  both  steel  and  graphite/epoxy  are  considered  in  the  simulations. 
The  plates  are  restrained  in  the  z-direction  at  the  outer  three  node  points  along 
the  upper  and  lower  surface  of  the  plates.  This  results  in  an  effective  plate 
diameter  of  267  mm.  The  finite  element  discretization  for  the  25.4  mm  thick 
plate  and  the  50.8  mm  thick  plate  are  shown  in  Figure  7.12  and  Figure  7.13, 


Figure  7.12  Finite-element  discretization  for  impingement  of  a 
spherical  shock  wave  on  a  25.4  mm  thick  plate. 

respectively.  The  origin  of  the  spherical  wave  is  at  R=0,  Z=-245  mm. 

Of  interest  in  these  simulations  is  both  the  early  time  response  as  well  as  the 
late  time  response.  Early  time  response,  is  taken  here  to  be  that  part  of  the 
response  which  is  governed  by  wave  propagation  in  the  through-thickness 
direction  of  the  plate.  The  through-thickness  stresses  are  responsible  for 
spallation  or  delamination  type  damage  in  composite  plates.  For  the  25.4  mm 
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Figure  7.13  Finite-element  discretization  for  impingement  of  a 
spherical  shock  wave  on  a  50.8  mm  thick  plate. 

this  early  time  is  on  the  order  of  50  -100  psec.  The  late  time-response  is  that 
portion  of  the  response  dominated  by  overall  flexure  of  the  plate.  It  is  this 
portion  of  the  response  that  results  in  matrix  cracking,  fiber  breakage  and 
overall  flexural  collapse  under  sufficient  loading.  Late-time  response  can  be 
estimated  by  considering  the  fundamental  frequency  of  the  plate  with  coupled 
fluid  mass.  For  the  configuration  considered  here  the  late-time  response 
occurs  in  the  range  of  0.1 -1.0  msec. 

For  the  early  time  response  of  the  25.4  mm  plate,  the  stress  history  at  two 
elements  within  the  plate  will  be  investigated.  The  first  of  these  is  at  R=1.71 
mm  and  Z  =0.668  mm,  which  essentially  gives  the  stress  transmitted  to  the 
plate  directly  below  the  origin  of  the  spherical  loading.  The  normalized  stress 
history  for  both  the  graphite/epoxy  plate  and  that  for  a  steel  plate  is  given  in 
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Figure  7.14.  As  seen  in  the  previous  subsection,  the  stress  wave  transmitted  to 


Figure  7.14  Normalized  stress  near  the  solid-fluid  interface  in  25.4 
mm  thick  graphite/epoxy  and  steel  plates  subjected  to  a  spherical  shock 
wave. 

the  steel  plate  is  significantly  higher  than  that  transmitted  to  the  graphite/ 
epoxy  plate.  An  integration  of  the  stress  history  over  time  is  a  measure  of  the 
momentum  density  at  that  location.  The  z-component  of  momenmm  density, 
Mj,  is  given  by 

=  (7.5) 

where  n,.  and  nj,  are  the  components  of  the  unit  normal  vector  in  the  radial  and 
through  thickness  direction,  respectively.  Considering  the  stress  state  at  the 
surface  of  the  plate,  the  shear  stress  is  zero  and  the  outer  normal  is  along  the 
z-direction  so  that  Equation  (7.5)  reduces  to  the  integration  of  the  normal 
stress  history  in  the  z-direction  over  time,  which  is  shown  in  Figure  7.14. 
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Performing  the  integration,  the  momentum  density  transferred  to  the  graphite/ 
epoxy  is  approximately  50%  of  that  transferred  to  the  steel.  The  reduced 
amount  of  momentum  density  transferred  to  the  plate  is  significant  since  the 
total  momentum  must  be  resisted  by  stresses  within  the  plate  during  the 
flexural  deformation. 


The  through  thickness  stress  (a^  is  also  shown  for  an  element  that  is  at  the 
plate  mid-thickness.  The  stress  history  for  both  the  steel  and  graphite/epoxy 
plates  are  shown  in  Figure  7.15.  Of  interest  here  is  that,  because  of  the  slower 


wave  speed  in  the  thickness  direction,  the  graphite/epoxy  plate  develops 
tensile  stresses  after  the  reflection  of  the  wave  from  the  rear  surface. 


The  late  time  response  is  manifested  by  the  flexural  deformation  of  the  plate. 
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high  local  stresses  in  the  region  of  the  restraint.  In  simulation  of  the  EBT  it 
would  be  preferable  to  use  a  contact  surface  support  over  the  length  of  plate 
that  is  actually  restrained.  The  radial  stress  history  at  the  rear  surface  of  the 
plate  is  shown  for  both  steel  and  graphite/epoxy  in  Figure  7.17.  The 


Figure  7.17  Normalized  stress  at  the  rear  surface  of  25.4  mm  thick 
graphite/epoxy  and  steel  plates  subjected  to  a  spherical  shock  wave. 

alternating  tensile  and  compressive  stresses  at  the  point  again  shows  that  the 
plate  is  vibrating  in  essentially  its  fundamental  mode.  The  maximum  stress 
level  reached  is  slightly  higher  in  the  steel.  Maximum  elastic  deflection  in  the 
graphitc/epoxy  plate  is  2.79  times  that  in  the  steel,  but  of  course,  the  weight  of 
the  steel  is  5  times  that  of  the  graphite/epoxy. 

The  response  of  the  50.8  mm  thick  graphite/epoxy  plate  is  considered  next. 
The  plate  has  the  same  boundary  conditions  and  is  subjected  to  the  same 
loading  as  the  25.4  mm  thick  plate.  The  normal  stress  at  the  plate  surface  is 
shown  in  Figure  7. 18.  Comparing  the  pulse  to  that  shown  in  Figure  7.13  for 
the  25.4  mm  graphite/epoxy  plate  shows  it  to  be  longer  due  to  the  longer  time 
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Figure  7.18  Normalized  stress  Oq  at  the  solid-fluid  interface  of  a  50.8 
mm  thick  graphite/epoxy  plate  subjected  to  a  spherical  shock  wave. 

for  the  wave  to  travel  to  the  rear  surface  and  back  to  the  front.  The  normal 
stress  at  the  mid-plane  of  the  plate  thickness  is  shown  in  Figure  7.19.  Of 
particular  importance  is  that  now  the  tensile  stresses  that  develop  are 
significantly  larger  than  those  that  develop  in  the  25.4  mm  thick  plate,  and  are 
approximately  0.75  times  the  magnitude  of  the  incoming  wave.  Thus  for 
thicknesses  larger  than  50.8  mm,  tensile  stresses  approaching  the  magnitude 
of  the  incoming  wave  would  be  expected.  The  late  time  response  of  the  50.8 
mm  plate  is  illustrated  through  the  sequence  of  deformed  geometry  plots 
shown  in  Figure  7.20.  The  maximum  deflection  in  the  50.8  mm  graphite/ 
epoxy  plate  is  0.358  times  that  of  the  25.4  mm  graphite/epoxy  plate  or 
approximately  the  same  as  that  of  the  25.4  mm  steel  plate.  The  radial  stress 
history  at  the  center  of  the  rear  surface  of  the  plate  is  shown  in  Figure  7.21. 
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Comparison  with  Figure  7.17  shows  the  maximum  stresses  to  be  one-half  that 
of  the  thinner  plate. _ 


*  0  30  60  90  120 
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Figure  7.19  Normalized  stress  Oq  at  the  midplane  of  a  50.8  mm  thick 
graphite/epoxy  plate  subjected  to  a  spherical  shock  wave. 
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Figure  7.21  Normalized  stress  Off  at  the  rear  surface  of  a  50.8 
graphite/epoxy  plate  subjected  to  spherical  shock  wave. 


41 


Vm.  DEVELOPMENT  OF  TTDE  DAMAGE  MODEL 


A  continuum  damage  model  (CDM)  is  developed  to  simulate  the  onset  of 
damage  and  the  subsequent  softening  effects  of  various  degradation 
mechanisms.  This  capability  will  be  useful  to  design  tests,  to  analyze  test  data, 
and  to  simulate  various  operating  or  threat  conditions.  There  is  a  fairly  large 
literature  in  CDM  and  some  of  that  work  is  followed  in  the  development  given 
here.  In  particular,  the  works  of  Davison  and  Stevens  [29],  Krajcinovic  [30], 
and  Talreja  [31],  with  appropriate  simplifications  and  extensions,  are  used  for 
guidance. 

A  distributed  or  continuous  damage  model  is  particularly  appropriate  for 
composite  materials  because  of  their  tendency  to  arrest  propagating  flaws, 
soften  and  redistribute  the  loads.  At  high  loading  rates,  even  delamination 
phenomena,  which  might  behave  as  single  fracture  planes  on  slow  loading, 
tend  to  disperse  into  a  distributed  pattern  of  small  cracks  which  cannot  run  as 
fast  as  the  loading  process  progresses  and  which  may  be  amenable  to 
modeling  by  a  CDM  approach.  The  model  described  here  is  linear  in  strain, 
but  nonlinear  in  damage  due  to  softening  effects. 

MQt£mU}£SQ£iBliQR 

This  damage  model  is  restricted  to  a  thick  laminated  composite  material,  a 
portion  of  which  is  shown  in  Figure  8.1,  with  no  fiber  reinforcement  in  the  1- 
direction  (through-thickness  and  generally  the  wave  propagation  direction) 
and  a  balanced  (isotropic)  arrangement  of  reinforcing  fibers  in  the  2,3-plane 
(in-plane).  The  material  is  assumed  to  be  exactly  transversely  isotropic  with 
respect  to  both  stiffnesses  and  developing  damage  with  the  2,3-plane  being  the 
plane  of  isotropy.  This  material  symmetry  is  attained  approximately  by  the 
usual  filament  wound  or  flat  layup  composites  with  balanced  0°,  ±60°  (as  in 
Figure  8.1)  or  0°,  90°,  ±45°  layups.  The  approximation  of  transverse  isotropy 
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is,  of  course,  much  better  for  stiffnesses,  but  is  also  imposed  on  damage 
development  in  this  treatment. 


Dnma^e  Description 

The  evolving  modes  of  damage  are  restricted  to  certain  orientations  of  brittle 
cracking  of  the  matrix  material.  A  simple  maximum  strain  criterion  is  used  for 
fiber  breakage  with  no  allowance  for  evolution.  The  matrix  damage  is 
characterized  by  the  vector  V=(V|,V2,V3)  referred  to  the  material  coordinate 
system  shown  in  Figure  8.1  with  the  indices  denoting  direction  of  the 
prevailing  normal  to  the  plane  of  the  cracks  and  the  magnitudes  denoting  a 
measure  of  the  severity  of  the  damage.  This  severity  might  take  the  form  of  a 
normalized  local  areal  density  of  crack  surfaces.  However,  in  the  present 
phenomenological  treatment,  damage  magnitude  is  interpreted  simply  in 
terms  of  fractional  reduction  of  certain  elastic  properties  of  the  material. 
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The  Vi  component  of  damage  is  intended  to  represent  a  network  of  matrix 
cracking  predominately  aligned  in  the  2,3-plane  (with  normals  in  the  1- 
direction)  which  do  not  cross  reinforcing  fibers  and  which  lead  ultimately  to 
delamination  failure.  It  is  intended  that  V|  damage  will  describe  spallation 
damage  due  to  through-thickness  waves  and  the  progression  of  damage 
leading  to  complete  separation  due  to  spallation.  A  state  of  Vj  damage  is 
shown  schematically  in  Figure  8.2. 


The  remaining  matrix  damage  is  characterized  by  a  combination  of  V2  and  V3 
damage.  This  is  due  to  matrix  cracks  which  may  traverse  or  lie  between 
reinforcing  fibers,  but  which  do  not  break  the  fibers.  This  damage  is  assumed 
to  occur  in  a  randomly  distributed  fashion  both  spatially  and  by  orientation  (by 
combinations  of  V2  and  V3)  of  normals  in  the  2,3-plane  such  that  the  degraded 
material  properties  are  left  transversely  isotropic.  Such  will  be  the  case  if  the 
average  degradation  is  independent  of  the  orientation  of  the  V2,  V3 
components  of  the  damage  vector.  This  is  assumed  to  be  the  case  and,  further, 
that  the  material  property  degradation  depends  only  on  the  magnitude 
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Vs  =  M  +  v|  (8.1) 

of  the  in-plane  components  of  the  damage  vector.  Thus,  the  in-plane  damage 
Vg  is  characterized  by  a  scalar  with  indefinite  direction,  except  that  it  is 
oriented  perpendicular  to  the  Vj  component  of  damage.  The  and  Vs 
damages  are  contrasted  also  with  regard  to  their  evolution  to  very  different  end 
states.  Whereas  the  V]  damage  evolves,  perhaps  catastrophically,  to  a 
delamination  or  spallation  with  complete  separation,  the  Vs  damage  evolves 
to  a  spatially  saturated  state  of  damage  defined  by  Vs=l  and  is  followed  by 
other  modes  of  damage  (either  fiber  breakage  or  further  delamination 
damage).  The  Vs  damage  is  shown  schematically  in  Figure  8.3. 


Material  Model 

Talreja’s  [31]  development  for  vector  damage  in  a  transversely  isotropic 
material  is  followed  here.  The  assumption  of  small  damage  is  relaxed  since  the 
full  evolution  of  damage  to  an  end  state  is  sought.  The  material  response  is 
further  restricted  to  axisymmetry  with  the  1-axis  being  the  axis  of  symmetry 
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and  with  the  2-  and  3-directions  now  being  identified  with  r  (radial)  and  0 
(circumferential)  cylindrical  coordinates.  Placing  small  strain  and  reflection 
symmetry  requirements  (independent  of  sign  of  Vj  and  Vs)  on  V,  leads  to  a 
Helmholtz  free  energy  function 

PV  =  f(Ij.l2....,l6)  (8-2) 

of  the  strain  and  damage  invariants 

”  ^11 

h  ~  ^22  ■*■^33 
I  = 


*3  ‘'22  ‘'33 

I  =  £^ 

M  “  ^12 

I  — 

I  —  +  V? 


(8.3) 


The  stress  tensor  components  result  from  the  derivatives  ajj=9(pv)/  9ejj  and  are 
expressed  in  the  usual  transversely  isotropic,  axisymmetric,  linear  constitutive 


form  as 


^11  ^11^12^12  ®  ^11 

^22  _  ^12^22^23  ®  ^22 

^33  ^12^23^22  ®  ^33 

0^2  0  0  0  £j2 


(8.4) 


where  the  elastic  coefficients  Cy  are  functions  of  the  damage  invariants 
l5=Vi^  and  facilitate  simple  assumptions  on  damage  dependence, 

(8.4)  is  inverted  to  the  following  compliance  form 


'22  _ 


"''12^^11 

-''12/^11 


-''12/^11 


^^23^^22 


-^12/^11 

-V23/E22 


(8.5) 


1/(2G,,)  o. 
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where 


(1  V23)Ejj 


1-Vj3-(2v22E,2)/E„ 

['-(vW)/Ei,1E;2 

(l+V23)[l-V23-(2vf2E22)/E„] 

^  _  _ ^12^22 _ 

I-V23-(2vW)/Ei1 

[V23  ~  (^12^22^  ^E] )  ]  E22 
(I+V23)  t^“''23“  (2v22E22)/Ejj] 


(8.6) 


C44  “  ^^12  • 

The  damage  dependence  is  imposed  on  the  compliances  in  (8.5)  since  they 
allow  a  more  intuitive  interpretation  of  damage  effect  than  do  the  elastic 
coefficients  Cy  in  (8.4).  This  can  be  assumed  in  a  general  functional  form,  to 
be  set  by  experimental  evidence.  However,  for  the  sake  of  demonstrating 
possible  damage  effects  on  constitutive  behavior,  the  following  simple 
damage  dependences  are  assumed: 

E„  =  (1-V?)E», 

E22  =  (•“®i'^|)E22 

G,2  =  (1-V^)  (l-a2V|)G®2  (8.7) 

v,2=  (I-v2)(I-a3V|)v«2 

V23  =  ( I  -  Vj)  V23 

where  the  superscript  “0”  denotes  virgin  properties  and  the  fractions  ()<ai<l, 
i=l,4  are  included  to  prevent  complete  loss  of  material  integrity  as  the 
saturation  state  Vg-^l  is  reached.  Integrity  may  remain  due  to  unbroken 
reinforcing  fibers  in  the  2,3-plane.  The  particular  functional  dependence  on  V  j 
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and  Vs  in  (8.7)  is  quite  arbitrary  except  for  the  fact  that  the  Eu  and  E22 
dependence  can  be  taken  to  define  Vj  and  Vs  and,  in  such  case,  the  remaining 
functions  should  be  determined  by  corresponding  experimental  evidence  of 
material  degradation  effects  on  those  other  material  properties.  TTie  elastic 
coefficient  dependence  on  damage  is  then  determined  by  combining  (8.7)  with 
(8.6). 

Dama2e  Evolution 

Rate  dependency  (actually,  time  dependency)  is  introduced  into  the  CDM 
process  by  way  of  damage  evolution  with  assumed  dependence  on  the  current 
state  of  damage,  some  overstress  above  a  current  threshold,  and  material 
properties  controlling  evolution  rates.  There  is  such  a  sparsity  of  data  on 
damage  evolution  that  guidance  depends  mostly  on  general  intuition  and 
awaits  definitive  experiments.  Curran,  et  al,  [32]  give  extensive  consideration 
to  evolution  (growth)  and  nucleation  of  voids  and  Krajcinovic  [30],  in  a 
review  article,  addresses  the  range  of  possible  kinetic  laws  for  damage 
evolution.  The  intended  application  here  is  to  composite  materials  which 
usually  abound  with  imperfection  sites  from  which  cracks  can  grow.  Thus, 
nucleation  of  new  cracks  is  not  addressed  here  and  the  focus  is  on  growth  or 
evolution  of  damage. 

Both  the  Vj  and  Vg  types  of  damage  are  assumed  to  be  governed  by  a 
threshold  of  the  form 

F(  g,  f  (V))  <  0  for  no  damage  growth  (8.8) 

>  0  for  damage  growth 

where  F  is  a  scalar  threshold  function,  g  is  the  current  stress  tensor,  f  is  an 
array  of  current  threshold  parameters,  which  is  a  function  of  V,  the  current 
damage  vector.  The  details  of  this  threshold  fimction  are  illustrated  with  the 
following  specific  example. 
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A  hyperbolic  threshold  surface  of  the  Mohr-Coulomb  type  is  taken  to  be 
dependent  on  the  stresses  o.  t  (tension  and  shear)  as  follows: 

F  (o,  f)  =  +  -  (f ^  -  a) /f^  (8.9) 

where  the  parameters  f  are  related  to  specific  growth  threshold  strengths  and 
the  Coulomb  friction  tangent  as 

“  ^2 

To  = 

These  tension  and  shear  growth  thresholds  and  the  Coulomb  friction  tangent 
in  (8.10)  are  assumed  to  be  functions  of  the  damage  Vj  (or  Vg).  The  threshold 
surface  F=0  from  (8.9)  is  shown  in  Figure  8.4  along  with  parameters  and  the 
shortest  distance  d  from  an  external  stress  state  (o,x)  to  the  threshold  surface. 

For  example  to  complete  the  threshold  specification,  the  Vj  delamination 
damage  threshold  is  postulated  to  depend  only  on  the  stress  components 


The  growth  threshold  stress  and  Coulomb  friction  tangent  in  (8.10)  are 
postulated  to  depend  on  the  damage  Vj  as 

^Gl  ~  “^I^^GIO 

“^Gl  "  (8-12) 

^G1  “  ^GIO'*'^!  “^GIO^ 

where  the  “0”  subscript  denotes  virgin  (Vi=0)  threshold  properties.  These 
simple  functions  are  designed  so  that  all  resistance  to  damage  is  lost  as  V^-^l 
and  so  that  the  friction  tangent  can  vary  linearly  with  as  damage 
progresses.  The  f  dependence  on  V  is  then  established  by  combining  (8.10) 
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Figure  8.4  Threshold  surface  for  the  onset  of  damage, 
and  (8.12).  I  he  evolution  (or  rate)  equations  tor  the  damage  are  required  to 

complete  this  CDM  constitutive  description.  For  the  Vj -damage,  it  is  further 
postulated  that  the  material  derivative  takes  the  form 

dVi 

=F,«i,,V,)  (8.13) 


where  dj  is  the  shortest  distance  in  the  o=Oii,  xi=oj2  stress  space  from  an 
exterior  stress  point  to  the  threshold  surface  F=0  shown  in  Figure  8.4.  Fj  in 
(8. 1 3)  should  vanish  for  di=0,  making  dV  j/dt=0  for  all  stress  points  interior  to 
or  on  the  threshold  surface.  For  purposes  of  illustration,  (8.13)  is  written  in  the 
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(8.14) 


specific  power  law  form 

dt  =  (“/‘’g.o) 
where  iij  is  a  positive  power  term  for  dj/oGiO’  dimensionless  stress 
distance,  tii  is  a  time  constant  governing  the  rate  magrnfjde,  and  the  1/(1-V|^) 
term  causes  an  acceleration  of  V  2  damage  to  complete  delamination  as  V  1 . 

In  general,  the  rate  of  damage  may  be  a  function  of  the  location  g  in  the  stress 
space  and  the  separation  from  the  threshold  surface.  For  example,  it  could  be 
expected  that  shear-induced  damage  under  compressive  stress  will  progress 
much  more  slowly  than  will  tensile  stress-induced  damage. 


The  Vs  in-plane  isotropic  damage  is  postulated  to  depend  only  on  the  stress 
components 


^  ~  2  ^^22'*’ ^33^ 


(0.1 5) 


-  U2 


^T2+4^''22-^33)  • 


These  quantities  are  derived  from  certain  invariants  of  the  stress  tensor  in 
absence  of  the  oi  1  component.  This  stress  component  has  been  ignored  for 
simplification  and  because  its  effect  is  a  sensitive  function  of  Poisson’s  ratio 
coupling  which  is  very  uncertain  at  this  level  of  approximation  of  material 
property  evolution.  Similarly,  the  C22  <^33  stress  components  were  ignored 

in  their  possible  effect  on  V|-damage  evolution  in  (8.11).  The  stress 


components  o,  x  in  (8.15)  are  used  in  the  threshold  condition  (8.9)  to  establish 
a  distance  ds  to  the  threshold  surface.  The  threshold  parameters  ?  are  again 
given  by  (8.10)  and  the  remaining  quantities,  to  set  the  Vg  evolution,  follow. 
Since  the  damage  consists  of  an  ever  denser  network  of  cracks  which  tends 
to  saturate,  the  evolution  equation  for  V5  is  postulated  to  be 


Fs(ds,  Vg) 


(8.16) 
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with  weak  dependence  on  Vs-  That  is,  no  acceleration  to  a  catastrophic  failure, 
as  in  the  case  of  Vi-^1.  In  fact,  for  purposes  of  illustration,  take 

3i  =  (VW 

with  no  explicit  dependence  on  V5  and  where  115  is  again  a  time  constant,  dg 
is  normalized  with  the  virgin  growth  threshold  stress  Ooso’  corresponding  to 
a  in  (8.15),  and  ng  is  a  positive  power  term  for  the  dimensionless  stress 
distance.  Saturation  can  be  forced  through  the  growth  thresholds  as 

^GS  "  ^GSO^  ^  ^ 

"gs  =  w  (8.18) 

‘(’gS  ”  'f’GSO‘'''(s(*t’GSl  “'(‘gSO^ 

SO  that  V5  damage  is  increasingly  difficult  to  drive  as  Vs->1.  This  is  only  one 
possible  form  of  the  threshold  functions.  In  the  following  section,  another 
threshold  form  for  Vg  damage  is  considered,  which  may  more  closely 
represent  the  available  experimental  data. 

The  third  mode  of  damage  involves  reinforcing  fiber  breakage  and  it  may  also, 
in  general,  be  rate  dependent  due  to  statistical  distributions  of  fiber  strengths. 
For  high  fiber  volume  fraction  and  load  aligned  with  the  fibers,  it  has  been 
found  that  there  is  very  little  rate  dependence  (see  Harding  and  Welsh  [13]). 
In  keeping  with  the  simplifications  and  assumptions  which  have  already  been 
made,  a  simple  maximum  strain  criterion  for  fiber  breakage  is  used  here.  It  is 
influenced  by  continuum  damage  evolution  only  through  the  softening  effects 
of  Vj  and  Vg  damage  upon  the  composite  deformation. 
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IX.  DAMAGE  MODEL  RESULTS  FOR  HOMOGENEOUS 
DEFORMATION 

The  damage  model,  as  developed  in  Section  Vin,  is  exercised  for  spatially 
homogeneous  strain  fields  and  for  various  strain  rates  and  transients. 
Predictions  are  compared  with  experimental  observations  to  set  damage 
material  properties  and  to  test  the  versatility  of  the  CDM  approach. 
Composites  of  interest  to  this  program  do  not  have  a  strong  coupling  between 
the  V I  and  modes  of  damage  due  to  weak  Poisson’s  ratio  coupling  between 
the  deformations  which  drive  each  of  the  damage  types.  As  a  result,  combined 
damage  is  not  expected  except  in  biaxial  or  triaxial  deformation  Helds  where 
each  mode  of  damage  is  individually  driven  by  appropriate  components  of  the 
deformation.  For  the  sake  of  simplicity,  and  to  utilize  available  experimental 
data,  some  simpler  uniaxial  deformation  fields  are  considered  in  this  section. 

The  first  case  to  be  considered  involves  various  constant  strain  rates  causing 
Vg  in-plane  damage.  Analyses  and  experiments  on  these  cases  are  useful  to 
help  isolate  the  effect  of  strain  rate  on  damage  and  deformation.  It  especially 
should  be  noted  that  the  current  CDM  has  no  explicit  dependence  on  rates. 
Rather,  the  rate  dependence  will  manifest  itself  through  time  required  for 
evolution  of  the  underlying  damaging  and  softening  processes. 

The  second  case  involves  V  j  damage.  Delamination,  or  V  j  damage,  is  difficult 
to  monitor  in  slow  rate  tests  and  is  usually  absent  or  catastrophic  such  that 
intermediate  states  are  rarely  seen.  Impact  and  spallation  tests  are  most 
appropriate  for  controlled  levels  of  delamination  damage.  These  tests  do 
involve  wide  variations  of  deformation  rates  and  may  be  subject  to  qualitative 
judgements  on  the  state  of  damage.  Simulations  of  deformation  histories 
expected  in  wave  traversals  are  used  to  generate  predicted  damage,  which  is 
then  compared  with  experimental  evidence  for  spallation  damage  thresholds. 
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Constant  Strain  Rate  Tests 

Recent  data  has  been  obtained  by  Behler,  Sikordki,  and  Staskewitsch  [14]  for 
graphite/epoxy  composites  subjected  to  a  wide  range  of  constant  strain  rate 
loadings.  These  experiments  were  done  in  a  cantilever  configuration  with 
impact  furnishing  the  higher  loading  rates  and  hydraulic  loading  for  the  lower 
rates.  This  arrangement  does  not  precisely  match  the  symmetries  required  for 
transversely  isotropic  Vg  damage  since  these  tests  are  approximately  uniaxial 
flexural  stress  loadings  with  matrix  cracking  expected  to  be  orientated 
predominately  with  normals  parallel  to  the  tensile  stress  direction.  Also,  these 
certainly  are  not  homogeneous  deformation  fields.  However,  the  comparison 
of  experiment  with  COM  prediction  is  made  anyway  because  it  is  thought  that 
the  general  behavior  of  damage  evolution  in  the  uniaxial  flexural  tests  may  be 
quite  similar  to  that  for  a  homogeneous  biaxial  deformation  field. 

The  material  of  most  interest  here,  which  was  tested  by  Behler,  et  al  [14],  is  a 
multidirectionally  reinforced  graphite/epoxy  consisting  of  a  balanced  0°,  90°, 
±45°  layup  with  T800  graphite  fibers  and  Vicotex  6376  epoxy.  The  test 
specimens  measured  6  mm  x  6  mm  x  45  mm.  The  ends  of  the  cantilever 
specimens  were  struck  with  a  hammer  device  on  a  flywheel  which  was 
hydraulically  driven  for  low  impact  velocities  (2x10'^  m/sec  to  1  m/sec)  and 
was  electrically  driven  for  higher  impact  velocities  (5  m/sec  to  50  m/sec).  The 
specimens  were  monitored  by  strain  gauges  on  the  compression  and  tension 
sides  and  resulting  strain  rates  ranged  from  5x10'^  sec*^  up  to  1200  sec'^ 
While  the  deformation  field  varies  through  the  thickness,  the  maximum  strain 
and  damage  is  expected  to  be  localized  about  the  outer  fibers  on  the  tension 
side  and  that  response  is  compared  with  the  CDM  predictions.  These 
experimental  results  are  reproduced  in  Figure  9.1.  Definite  strengthening  and 
stiffening  with  increased  loading  rate  is  seen,  however,  little  variation  in 
ultimate  strain  is  indicated.  The  predominate  softening  mechanism  is  assumed 
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Figure  9.1  Flexural  impact  data  for  graphite/epoxy  (Behler,  et  al[14]). 
to  be  due  to  matnx  cracking  with  only  the  linal  curvature  near  the  breaking 
point  being  due  to  some  evolutionary  behavior  in  the  fiber  fracture. 


The  damage  model  in  Section  Vin  is  applied  to  biaxial  deformation  with 
011=0  and  £22=^33  *^®se  strains  being  at  the  constant  rates  shown  in 

Figure  9.1 .  The  first  attempt  is  with  a  model  that  has  a  finite  threshold  for  the 
virgin  state  of  no  damage  (Vs=0)  as  in  (8.18)  with  oqso  being  the  initial 
threshold  and,  beyond  this  threshold,  the  damage  evolves  according  to 
(8.17).The  pertinent  material  properties  are  given  in  Table  9.1. 


Table  9.1  Material  Properties  for  Damage 
Threshold  Model  A 

E?,  (MPa) 

10,000 

E?,  (MPa) 

75,000 

vO 

0.04 

vO 

0.30 

ai=a3 

0.75 

«4 

0.95 

"S 

1 

hs  (sec) 

10*^ 

^Gso  (MPa) 

950 

The  damage  Threshold  Model  A  refere  to  the  threshold  stress  in  the  in  plane 
direction  given  by  the  function  of  damage  Vg  as  in  (8.18)  or 

These  properties  and  this  threshold  model  are  combined  with  the  evolution 
equation  (8.17),  where  for  the  assumed  symmetries  in  this  deformation, 
dg  =  Maximum  [0,  <5^2  - 

The  elastic  properties  in  Table  9.1(E®j ,  £22*  v®2.  and  V23)  are  chosen  to  be 
representative  of  a  balanced  layup  of  graphite/epoxy.  The  fractional 
reductions  (C],  03,  and  04)  are  chosen  very  arbitrarily  to  attempt  to  reflect  the 
elastic  integrity  remaining  in  a  saturated  matrix  cracked  state  due  to  the  intact 
fiber  reinforcements.  The  exponent  power  was  left  at  ns=l  under  the 
assumption  that  the  stress,  over  a  threshold,  drives  the  damage.  This  left  ris 
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ogso  ^  the  adjustable  parameters  in  the  CDM  to  attempt  to  match  the 
experimental  results  in  Figure  9.1.  Ihe  approach  taken  was  to  adjust  the  virgin 
threshold  to  account  for  the  apparent  breakpoint  for  softening  in  the  fastest 
three  loading  paths.  That  is,  the  950  MPa  value  in  Table  9.1  represents  the 
022=^33  stress  point,  on  initial  loading,  below  which  absolutely  no  Vg  damage 
occurs.  Then  the  rate  parameter  iig  was  adjusted  to  best  fit  the  softening 
behavior. 

The  results  are  shown  in  Figure  9.2  where  the  CDM  prediction  with  Threshold 
Model  A  (defined  by  Table  9.1  and  Equation  (9.1))  are  compared  individually 
with  each  of  the  stress-strain  curves  from  Figure  9.1.  This  comparison  is  quite 
favorable  at  high  strain  rates  with  the  threshold  and  softening  behavior 
allowing  the  predicted  stress-strain  curve  to  somewhat  follow  the 
experimental  curves.  Especially  the  1200  sec’^  comparison  shows  a  small 
amount  of  softening  at  the  upper  end,  which  follows  the  experimental  data 
nicely.  These  calculations  were  run  out  to  a  fixed  strain  of  0.015  (for  a  rate 
independent  fiber  failure  criterion)  and  terminated.  The  final  stages  of 
softening  on  the  12(X)  sec'^  curve  may  be  predominately  due  to  evolution  of 
fiber  breakage,  which  is  not  accounted  for  in  this  CDM.  At  rates  below  0.85 
sec'\  the  predicted  curves  become  almost  exactly  bilinear  and  invariant  with 
respect  to  rate.  At  these  slow  rates,  the  threshold  moves  with  the  loading  and 
the  damage  evolution  is  in  a  condition  of  quasi-static  equilibrium.  That  is,  the 
damage  will  cease  to  increase  immediately  if  the  loading  ceases  to  increase. 
The  nonzero  virgin  threshold  stress  is  seen  to  fail  rather  badly  for  the  lowest 
(0.005  sec'^)  rate  of  loading,  which  exhibits  apparent  softening  for  any 
nonzero  loads.  These  observations  form  the  basis  for  another  threshold  model, 
which  follows,  and  which  attempts  to  more  directly  fit  the  raw  data  from 
Figure  9.1. 
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Figure  9.2  Threshold  Model  A  correlations  with  data  of  Figure  9.1. 


Since  it  is  expected  that  the  loading  stress  and  the  threshold  stress 
approximately  coincide  for  low  loading  rates,  the  0.005  sec"^  data  of  Figure 
9.1  was  fit  directly  by  a  least  squares  routine  to  the  threshold  stress.  This  was 
done  by  estimating  the  actual  damage  from  the  reduced  modulus  on  this 
loading  path  and  then  expressing  the  stress  as  a  function  of  damage  (rather 
than  strain).  The  required  material  properties  to  define  this  model  are  listed  in 
Table  9.2 


Table  9.2 

Material  properties  for  Damage 
Threshold  Model  B 

E®,  (MPa) 

10,000 

E?,  (MPa) 

75,000 

V® 

0.04 

V® 

^71  . 

0.30 

0.75 

«4 

0.95 

Os 

1 

hs  (sec) 

0.0014 

ogso  (MPa) 

170 

This  Threshold  Model  B  has  an  in-plane  threshold  function  given  by 

This  form  does  not  saturate  as  Vj-^l  (where  ogs"^3ogso)>  ^^es 

closely  follow  022  for  the  slowest  (0.005  sec'*)  rate  loading,  ogso  adjusted 
to  achieve  this  agreement  and  then  the  rate  constant  115  was  adjusted  to  best 
match  the  56  sec'^  test  data.  In  Figure  9.3,  comparisons  with  experiment  are 
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shown. 


This  Threshold  Model  B  is  seen  in  to  fit  a  wide  range  of  data  very  well. 
Deficiencies  are  present  at  the  later  stages  of  loading  where  fiber  failure  is 
probably  dominating  the  detailed  softening  and,  as  in  Figure  9.2,  for  all  strain 
rates  at  least  below  0.85  sec’^  the  model  predictions  are  invariant  with  respect 
to  rate.  As  a  result,  the  agreement  is  quite  poor  for  the  0.85  sec’^  data.  Possible 
explanations  are  1)  that  another  slow  mechanism  of  softening  should  be 
present  as  an  additive  factor  in  oqs  (perhaps  viscoelastic  deformation)  or  2) 
that  the  data  itself  should  be  rate  invariant  at  these  low  rates.  Conclusions  will 
have  to  await  further  experimentation  and  study. 

Wave  and  Spallation  Simulation 

The  most  appropriate  experiments  for  validation  of  this  damage  model  for 
delamination  or  spallation  damage,  of  the  type  shown  in  Figure  8.2,  are  planar 
impact  by  a  flyer  plate  onto  a  material  sample  with  no  backing.  This 
configuration  gives  rise  to  a  transient  tensile  pulse  due  to  the  release  wave 
emanating  from  the  stress  free  back  face.  Such  tests  and  their  assessments  are 
in  progress  under  this  program  and  will  be  reported  at  a  later  time.  For  the 
purpose  of  an  immediate  comparison  of  spallation  predictions  with 
experiments,  Roylance  [3]  gave  qualitative  results  of  spallation  thresholds  in 
graphite/epoxy  for  variation  of  incident  pulse  shapes  and  those  conclusions 
will  be  tested  here  against  simulated  wave  deformations  in  the  form  of 
homogeneous  fields.  Roylance ’s  results  are  not  detailed  enough  to  allow 
determination  of  damage  variables,  as  was  done  in  the  preceding  subsection 
for  the  V5  damage,  and  so  only  a  cursory  application  of  the  model  with 
estimated  properties  will  be  given  here  to  demonstrate  the  behavior  of 
spallation  or  Vj  damage  predictions  under  transient  loading. 

Roylance  in  [3]  considered  a  variety  of  graphite  reinforcing  fibers  (Pluton, 
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Hercules,  and  Celanese)  in  an  American  Cyanamid  BP907  epoxy  and  reported 
mainly  longitudinal  (fiber  direction)  properties.  Transverse  direction 
properties  (1 -direction  from  Figure  8.1)  are  of  the  most  interest  here  for 
spallation  damage  and,  hence,  the  epoxy  bonding  properties,  rather  than  the 
fiber  properties  are  thought  to  be  most  important.  Indeed,  Roylance  found  no 
dependence  of  spallation  threshold  on  the  graphite  fiber  type,  even  though  the 
longitudinal  fiber  strengths  and  moduli  were  reported  in  [3]  to  vary 
considerably.  As  a  result,  the  choice  was  made  here  to  use  representative 
graphite/epoxy  properties  with  the  realization  that  the  following  results  are 
most  sensitive  to  the  transverse  directional  properties.  These  properties  are 
listed  in  Table  9.3  along  with  Vj  damage  properties.  These  damage  properties 


Table  9.3  Material  properties  for  spallation 
simulation 

E?,  (MPa) 

10,300 

(MPa) 

69,000 

vO 

0.04 

vO 

0.30 

01=03 

0.75 

O4 

0.95 

ni 

2 

hi  (sec) 

4x10-^ 

^Gio  (MPa) 

50 

are  being  found  to  work  well  in  correlating  impact  test  data  in  ongoing  and  yet 
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to  be  reported  work  and  they  will  be  used  here  also  for  the  purpose  of 
illustrating  the  model  and  comparing  with  Roylance’s  results,  to  the  extent  of 
interpreting  spallation  thresholds. 

Symmetries  in  the  configuration  of  interest  allow  a  reduction  to  uniaxial  strain 
in  the  1 -direction  with  the  first  of  (8.12)  specifying  the  tensile  threshold  oq! 
and  (8. 14)  the  evolution  equation  for  V  j  damage.  In  fact,  the  symmetries  allow 
the  distance  from  the  damage  threshold  to  be  written  as  di=:Max[0,aii-aG2]. 
The  homogeneous  deformations  are  specified  through  the  only  nonzero  strain 
component  cn  as  a  function  of  time  for  all  points  in  space. 

The  first  case  considered  is  various  constant  strain  rates  from  a  state  of  zero 
deformation  into  tension.  This  is  done  to  illustrate  the  rate-dependent  behavior 
of  the  damage  model  and  to  compare  with  Figures  9.2  and  9.3  for  Vs  damage. 
These  results  are  seen  in  Figure  9.4  where  both  the  stress  versus  strain  and 
damage  versus  strain  are  seen  for  a  wide  range  of  constant  strain  rates. 
Obviously  both  tiie  stress  and  strain  at  failure  are  greatly  increased  as  the 
strain  rate  is  increased.  The  interpretation,  in  the  context  of  this  model,  is  that, 
rather  than  a  pure  rate  dependence,  damage  evolution  is  lagging  the 
application  of  load  for  the  more  rapid  deformation.  It  is  also  observed  that 
ultimate  failure  is  identified  as  complete  loss  of  load  carrying  capability 
through  loss  of  modulus  Eu  in  (8.7)  as  Vj-^l.  At  low  strain  rates  (10  and  100 
sec'^  in  Figure  9.4),  this  model  is  seen  to  be  a  simplistic  ultimate  strength 
criterion  with  the  load  carrying  capacity  falling  rapidly  for  02i>ogi. 

The  test  conditions  in  Roy  lance’s  paper  [3]  are  simulated  with  homogeneous 
deformation  fields  by  constructing  rectangular  and  triangular  strain  pulses. 
The  compressive  incident  strain  pulses  were  then  augmented  with  the  tensile 
release  pulse  from  a  presumed  free  back  face  and  superposed  as  if  on  the 
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Strain  Rate  10  lO^ 


Figure  9.4  Vj  damage  at  various  constant  strain  rates. 

spallation  plane,  defined  as  that  position  nearest  the  back  face  which  first  feels 
maximum  tension.  This  approach  is  somewhat  simplistic  and  it  drives  the 
damage  without  feedback  to  the  deformation  field,  but  it  does  provide  a 
prediction  for  softening  of  the  elastic  moduli  and  it  is  much  easier  to 
implement  than  is  a  full-blown  wave  propagation  and  damage  analysis.  Using 
the  constructed  strain  pulses,  the  stress  and  damage  histories  are  compute  1 
with  the  model.  Wave  propagation  calculations  are  not  really  justified  or 
appropriate  for  comparison  with  experimental  results  unless  fairly  extensive 
instrumentation  and  data  acquisition  has  provided  quantitative  rather  than 
qualitative  data.  Some  of  these  more  extensive  analyses  are  presented  in 
Section  X. 

Ro' lance  considered  exploding  foil  experiments  on  graphite/epoxies  with 
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pulse  shapes  being  approximately  triangular  if  the  pressure  pulse  from  the 
explosion  directly  encountered  the  material  and  rectangular  if  the  pressure 
pulse  drove  a  flyer  plate  into  the  sample  material.  The  three  pulses  associated 
with  spallation  thresholds  were  rectangular  of  amplitude  90  MPa  and  duration 
1.08  psec,  rectangular  of  amplitude  125  MPa  and  duration  0.67  psec,  arid 
triangular  with  amplitude  143  MPa  and  duration  1.42  psec.  These  pulses  all 
deliver  an  impulse  to  the  material  of  approximately  1  kilotap  (1  kilclap=l 
kilobar*psec=100  Pa«psec).  Specifically,  if  the  pulse  rises  (and  decays  for  the 
rectangular  shapes)  are  instantaneous,  which  will  be  assumed  here,  then  the 
impulses  are  0.972, 0.837,  and  1.015  kilotaps,  respectively,  in  the  order  listed 
above.  These  stress  pulses  were  converted  to  strain  deformation  pulses 
assuming  the  virgin  elastic  properties  given  in  Table  9.3  and  the  uniaxial  strain 
symmetry  condition.  The  strain  pulses  arc  shown  in  Figure  9.5  with  the  tensile 
release  wave  appended  as  expected  for  the  deformation  on  the  spallation 
plane.  Also  shown  in  Figure  9.5  are  the  corresponding  predicted  stress  pulses 
with  a  slight  bit  of  stress  relaxation  due  to  damage  softening  during  the  tensile 
portion  of  the  deformation.  The  Vj  damage  as  a  function  ot  time  is  shown  as 
the  third  transient  record  in  Figure  9.5.  Since  the  initial  damage  threshold 
ogio=50  MPa  from  Table  9.3  is  obviously  exceeded  in  each  of  the 
deformations,  it  is  not  surprising  that  damage  does  occur.  It  is  observed  that 
the  predicted  reduction  in  the  Eu  modulus,  from  Equation  (8.7),  is  V|^  and 
only  ranges  from  approximately  2-12%.  This  is,  perhaps,  in  the  appropriate 
range  to  be  identified  as  spallation  threshold,  as  was  done  by  Roy  lance  in  [3]. 
Roylance  had  speculated  that  the  approximately  1  kilotap  impulse  delivered  in 
each  case  was  indicative  of  an  impulse  criterion  for  spallation  threshold.  If  the 
present  damage  model  is  representative  of  real  material  damage  behavior,  then 
impulse  is  not  an  appropriate  criterion  (accounting  for  the  fairly  wide 
variations  in  damage  at  approximately  equal  impulse)  and,  instead,  the  load 
sustained  by  the  material  above  damage  threshold  controls  the  evolution  of 
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Figure  9.5  Strain,  stress,  and  damage  histories  predicted  for 
approximately  1  kilotap  pulses  of  various  shapes. 

damage.  Conclusions  on  these  questions  await  further  test  results  and 
correlations  with  damage  models.  It  is  concluded  that  softening  is  a  very 
important  part  of  the  total  dynamic  material  response  and  that  the  CDM 
approach  is  a  convenient  and  appropriate  way  to  handle  these  situations. 

For  the  purpose  of  illustrating  more  substantial  damage  in  a  dynamic 
deformation,  the  strain  deformation  pulses  of  Figure  9.5  were  magnified  and 
the  softening  and  damage  were  ^calculated.  These  results  are  shown  in  Figure 
9.6  for  an  impulse  delivered  of  approximately  twice  that  of  Figure  9.5.  The 
stress  magnitudes  were  adjusted  such  that  the  compressive  part  of  the  pulse  for 
all  three  cases  delivers  exactly  2  kilotaps.  This  is  done  by  making  the  1 .08  psec 
duration  rectangular  pulse  have  magnitude  185.2  MPa,  the  0,67  psec  pulse 
have  magnitude  298.5  MPa,  and  the  1 .42  psec  triangular  pulse  have  magnitude 
28 1 .7  MPa.  It  is  seen  in  Figure  9.6  that  now  the  softening  and  damage  are  very 
dramatic  with  rapid  softening  immediately  upon  application  of  substantial 
tensile  overstress  and  damage  levels  almost  completely  to  Vj=l  where  the 
material  has  lost  all  stiffness.  This  inability  to  support  large  tensile  stresses  on 
the  spallation  plane  is  what  alters  the  expected  (no  loss  of  integrity)  back-face 
wave  response  (usually  in  the  form  of  back  face  particle  velocity)  and  forms 
what  is  called  a  spallation  signature. 
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Figure  9.6  Strain,  stress,  and  damage  histories  predicted  for  2  kilotap 
pulses  of  various  shapes. 


X.  PRELIMINARY  APPLICATION  OF  DAMAGE  MODEL  TO  TWO 
DIMENSIONAL  PROBLEMS 

In  order  to  illustrate  the  application  of  the  damage  model  to  two  dimensional 
axisymmetric  problems,  the  response  of  the  graphite/epoxy  plate  with 
diameter  267  mm  subjected  to  shock  loading  from  underwater  explosive  will 
be  reconsidered.  The  material  propeities  to  be  used  in  the  simulations  will  be 
those  given  in  Section  IV  for  the  virgin  material  along  with  the  damage 
parameters  used  in  Threshold  Model  B  of  Section  DC.  In  addition  to  the 
parameters  considered  there,  the  fully  two-dimensional  treatment  requires 
specification  of  several  additional  parameters  related  to  damage  under  shear 
stresses.  The  parameters  used  here  are  estimated  and  intended  principally  to 
demonstrate  application  of  the  material  model.  A  complete  summary  of  the 
parameters  used  is  given  in  Table  10.1. 

As  discussed  in  Section  VI,  the  use  of  local  zero  displacement  boundary 
conditions  along  the  outer  diameter  of  the  plate  results  in  very  large  stresses 
concentrated  in  that  area.  In  order  to  avoid  having  damage  occur  only  along 
the  outer  edge,  the  plate  will  be  taken  to  be  free  in  the  z-direction.  This  results 
in  a  somewhat  higher  fundamental  frequency  of  the  plate  than  the  simply 
supported  condition  considered  in  Section  VI  and  lower  radial  stress  for  a 
given  loading  history.  Thus,  somewhat  higher  pressures  are  required  to  cause 
material  damage  than  would  be  needed  in  a  plate  supported  along  the  outer 
diameter. 

The  introduction  of  damage  into  the  material  constitutive  behavior  causes 
nonlinearities  in  the  response  so  that  use  of  nondimensionalized  parameters  in 
examining  results  as  was  done  in  Section  VI  is  not  possible.  The  pressure 
history  used  here  is  one  selected  to  result  in  a  level  of  damage  in  the  plate  to 
illustrate  the  application  of  the  model.  The  pressure  distribution  of  the  incident 
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Table  10.1  Material  Properties  for  Two 
Dimensional  Problems 


E®,  (MPa) 

10,300 

(MPa) 

69,000 

vO 

0.04 

vO 

0.30 

G?,  (MPa) 

7,170 

ai=a2=«3 

0.75 

04 

0.95 

<^GSO  (MPa) 

170 

tosoCMPa) 

340 

"s 

1 

71s  (sec) 

0.0014 

<t>GSl=<l>GSO 

1.0 

<^Gio  (MPa) 

50 

''Gio  (MPa) 

50 

"1 

2 

Til  (sec) 

4x10*^ 

4>gii=<I>gio 

0.5 

^ult 

0.015 

I 


70 


wave  applied  at  the  surface  of  the  fluid  is  chosen  to  produce  an  incident 
pressure  at  the  surface  of  the  plate  at  R  =  0  of  160  MPa  with  a  characteristic 
time  0|  of  0.1 130  msec.  The  radial  distribution  of  P  and  6  is  then  determined 
by  Equations  (5.1)  and  (5.3)  respectively.  In  the  first  simulation  the  50.8  mm 
thick  plate  is  considered.  Because  of  the  presence  of  nonlinearity,  results  are 
presented  in  terms  of  the  z-direction  displacement  history  rather  than  stress 
history.  In  order  to  eliminate  the  rigid  body  motion  component  of 
displacement  associated  with  the  free  plate,  the  quantity  that  is  presented  is  the 
displacement  of  the  center  of  the  rear  surface  minus  the  displacement  at  the 
edge  of  the  rear  surface.  The  effect  of  damage  is  seen  by  a  comparison  to  the 
plate  response  with  the  absence  of  damage.  Figure  10.1  shows  the  computed 


TIME  no"” 


Figure  10.1  Displacement  history  (in  meters)  in  50.8  mm  thick 
graphite/epoxy  plate  subjected  to  spherical  shock  wave  (with  damage 
— ,  without  damage . ). 
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displacement  history  with  and  without  the  effects  of  damage.  As  seen  in  the 
figure,  the  response  of  the  plate  for  the  two  computations  is  almost  identical, 
indicating  minimal  amounts  of  damage.  The  damage  contours  themselves  are 
shown  in  Figures  10.2  and  10.3  at  a  time  of  500  ps.  The  Vi-damage  in  Figure 


10.2  is  concentrated  near  the  z-axis,  where  the  direct  through  thickness  tensile 
stresses  are  greatest  and  shear  stress  is  minimal  indicating  that  the 
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delamination  damage  occurs  from  the  stress  waves  propagating  through  the 
thickness  of  the  material.  The  contour  shown  is  for  a  value  of  0.25,  which 
represents  a  sub-spallation  condition,  with  some  reduction  in  the  material 
stiffness.  The  Vg-damage  contour  is  shown  in  Figure  10.3.  As  expected  the 
damage  occurs  at  the  upper  and  lower  surface  of  the  plate  near  the  center  as 
the  plate  oscillates,  producing  alternating  tension  and  compression  on  each 
face.  Again  though,  the  level  of  the  damage  is  fairly  low  so  that  the  reduction 
of  stiffness  is  not  substantial,  which  is  consistent  with  the  minor  differences 
noted  in  the  displacement  history. 

In  the  next  simulation  the  effect  of  a  higher  level  of  damage  is  considered.  A 
25.4  mm  thick  plate  subjected  to  a  shock  wave  of  twice  the  amplitude  as  the 
previous  simulation,  but  retaining  the  same  spatial  and  temporal  distribution, 
will  be  investigated.  The  same  measure  of  displacement  will  be  used  to  study 
the  effect  of  damage  on  the  plate  response.  The  displacement  history  for  the 
simulation  is  shown  in  Figure  10.4.  Here  it  can  be  seen  that  the  damage  has 
greatly  influenced  the  response  of  the  plate  after  the  initial  oscillation.  The 
stiffness  has  been  reduced  to  such  extent  that  the  displacements  at  the  time  of 
400  psec  for  the  plate  with  and  without  damage  are  almost  completely  out  of 
phase. 

The  distribution  of  damage  is  examined  at  both  early  and  late  time.  The  Vj- 
damage  distribution  is  shown  in  Figure  10.5  at  a  time  of  125  psec.  Contours 


Figure  10,5  Vj  damage  contours  at  125  psec  (N=0.95, 0=0.25). 
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Figure  10.4  Displacement  history  (in  meters)  in  25.4  mm  thick 
graphite/epoxy  plate  subjected  to  spherical  shock  wave  (with  damage 
— f  without  damage ). 

csz  shown  for  two  levels  of  damage,  one  representing  onset  or  partial  damage 


and  one  representing  delamination.  As  seen  in  the  figure,  there  is  no  damage 


along  the  z-axis,  indicating  that  the  damage  is  a  result  of  shear  stresses  not  the 


through-thickness  tensile  stresses  as  was  the  case  in  the  50.8  mm  thick  plate. 


The  contours  of  the  V^-damage  are  shown  in  Figure  10.6.  At  this  time  tension 
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has  only  occurred  on  the  rear  surface  of  the  plate,  thus  there  is  no  Vg-damage 
on  the  front  surface.  The  Vj  and  Vg  damage  contours  at  500  psec  are  shown 
in  Figures  10.7  and  10.8,  respectively.  Examination  of  the  damage  contours 


shows  that  the  delamination  or  Vj  damage  extends  almost  entirely  across  the 
radius  of  the  plate  (in  effect  saying  that  complete  delammation  has  occurred), 
thus  preventing  shear  transfer  through  the  thickness. 

The  final  example  to  be  considered  is  for  the  50.8  mm  thick  plate  subjected  to 
320  MPa  with  the  same  0,=O.l  130  msec.  The  Vj -damage  contours  are  shown 
at  the  time  of  500  psec  in  Figure  10.9.  As  seen  in  Figure  10.9,  there  is  a  large 
region  of  delamination  damage  near  the  z  axis  and  also  one  approximately  at 
two  thirds  of  the  plate  radius.  The  first  of  these  is  attributed  to  the  through- 
thickness  tensile  stresses  that  develop  while  the  second  region  is  attributed  to 
shear  stresses  since  this  is  a  region  of  inflection  in  the  plate  bending.  This  can 
also  be  seen  by  examining  the  deformed  geometry  of  the  plate  in  Figure  10.10 
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at  the  time  of  500  nsec.  The  central  portion  of  the  plate  shows  a  region 
disturbed  by  damage  corresponding  to  the  first  contour  of  critical  damage.  As 
Vi->1  several  of  the  material  coefficients  go  to  zero,  which  permits  the 
material  to  deform  extensively  under  very  small  loads,  resulting  in  the 
distorted  region  seen  in  the  figure.  The  other  region  of  high  V  j  damage  is 
manifested  by  the  large  shear  distortion  seen  across  the  thickness  of  the  plate. 
The  softening  of  the  material  with  damage  results  in  a  localized  band  in  which 
the  distortion  takes  place  somewhat  analogous  to  the  shear  bands  seen  in 
metals  under  shear  deformation. 
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XL  SUMMARY 


This  report  has  concentrated  on  the  early  time  material  response  of  thick 
composite  plates  by  presenting  the  development  and  demonstration  of 
techniques  and  models  appropriate  for  response  predictions  and  experimental 
correlation.  Various  levels  of  theoretical  and  numerical  prediction  techniques 
are  compared  with  available  experimental  data  and  with  each  other.  The  time 
and  space  regimes  of  interest  for  shock  wave  damage  clearly  include  both 
body  wave  response  and  early  flexural  and  shear  vibration  modes  when  two- 
dimensional  plate  configurations  are  considered.  Finite-element  calculations 
are  presented  which  show  direct  through-thickness  wave  arrivals  which 
transition  into  flexural  vibrations  as  time  progresses.  Specific  findings  are 
listed  below. 


The  results  presented  indicate  that  the  use  of  continuum  finite 
elements  in  conjunction  with  silent  boundary  conditions  to 
represent  the  fluid  is  an  appropriate  modelling  technique,  at  least 
for  two-dimensional  problems.  The  technique  allows  for  the 
effects  of  material  nonlinearity  and  cavitation  to  be  included 
without  a  tremendous  computational  penalty.  The  results  indicate, 
however,  that  a  more  refined  mesh  than  that  considered  here 
should  be  used  for  comparing  computational  results  to  actual 
experiments  or  for  predictions  of  component  response,  particularly 
where  through  thickness  stress  is  important.  For  predominately 
flexural  behavior  the  modelling  technique  employed  gave  results 
which  agreed  extremely  well  with  analytical  results. 

The  use  of  zero  displacement  boundary  conditions  at  the  edge  of 
the  plates  results  in  very  high  local  stresses.  For  cases  of  modelling 
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a  composite  material  EBT,  an  attempt  to  simulate  the  actual  plate 
support  conditions,  as  much  as  is  feasible,  should  be  made.  Use  of 
contact  surfaces  may  be  beneficial. 

Results  of  calculations  of  spherical  shock  waves  on  steel  and 
graphite/epoxy  plates  show  that  he  momentum  density 
transmitted  to  the  graphite  epoxy  plate  is  approximately  50%  less 
than  that  transmitted  to  the  steel  plate  due  to  its  lower  mechanical 
impedance.  This  is  significant  since  this  momentum  density  must 
be  resisted  by  the  radial  stresses  in  the  plate.  The  calculations  also 
show  that  for  plates  of  267  mm  diameter  and  25.4  mm  thickness 
the  displacement  of  the  graphite/epoxy  is  approximately  twice  that 
of  the  steel,  but  at  one  fifth  the  weight,  which  demonstrates  a 
distinct  advantage  of  the  composite  material.  However,  due  to  the 
relatively  slow  wave  speed  in  the  through  thickness  direction  of  the 
composite  plate,  tensile  stresses  develop  even  under  spherical 
waves.  These  tensile  stresses  are  particular  significant  due  to  the 
low  transverse  strength  of  the  composite,  indicating  that  spallation 
in  the  thick  composite  is  a  possible  damage  mode. 

A  continuum  damage  model  was  developed  for  application  to  thick 
composite  materials  subjected  to  high-rate  dynamic  loading.  The 
model  is  capable  of  describing  the  rate  dependent,  nonlinear 
behavior  of  the  composite  material.  The  model  compares 
favorably  to  the  limited  available  experimental  data  at  high  rates 
for  cases  of  idealized  hoinogeneous  deformation. 

Preliminary  application  of  the  model  to  multi-dimensional 
problems  was  performed  to  illustrate  its  capability.  The  results 
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obtained  appear  reasonable,  but  were  not  compared  to  experiments 
due  to  lack  of  available  data.  There  also  are  a  number  of 
parameters  in  the  model  that  are  difficult  to  determine  and  thus  far 
have  been  estimated,  with  some  parametric  studies  conducted  to 
assess  tlie  sensitivity  of  the  results.  Even  with  these  difficulties  the 
model  appears  to  be  a  very  promising  step  to  describe  the  material 
response  under  these  transient  loadings. 


This  work  is  supported  by  the  DARPA  Naval  Technology  Office  and  that 
support  is  greatly  appreciated. 
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